1 Setup

1.1 Load packages

library("knitr")       # for knitting stuff
library("kableExtra")  # for markdown tables
library("lme4")        # for linear mixed effects models
library("broom")       # for tidying model results 
library("brms")        # Bayesian regression with Stan
library("corrr")       # for tidy correlation matrix
library("xtable")      # for latex tables
library("jsonlite")    # for reading json files
library("janitor")     # cleans stuff
library("Hmisc")       # bootstrapped confidence intervals
library("tidybayes")   # for Bayesian data analysis
library("png")         # adding pngs to images
library("grid")        # functions for dealing with images 
library("plotly")      # 3D scatter plot 
library("egg")         # for geom_custom
# library("ggtern")      # for ternary plots 
library("patchwork")   # for figure panels
library("tidyverse")   # for wrangling, plotting, etc. 

1.2 Helper functions

# setting some knitr options
opts_chunk$set(comment = "",
               fig.show = "hold")

# set the ggplot theme
theme_set(theme_classic())

# function for printing out html or latex tables 
print_table = function(data, format = "html", digits = 2){
  if(format == "html"){
    data %>% 
      kable(digits = digits) %>% 
      kable_styling()
  }else if(format == "latex"){
    data %>% 
      xtable(digits = digits) %>%
      print(include.rownames = F,
            booktabs = T)
  }
}

# root mean squared error
rmse = function(x, y){
  return(sqrt(mean((x - y)^2)))
}

actual_counterfactual_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = clipindex,
                            y = value,
                            fill = colorindex,
                            group = model)) +
    stat_summary(geom = "bar",
                 fun = "mean",
                 color = "black",
                 position = position_dodge(0.7),
                 width = 0.7) +
    stat_summary(geom = "linerange",
                 fun.data = "mean_cl_boot",
                 position = position_dodge(0.7),
                 size = 1) +
    geom_point(data = x %>% mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(dodge.width = 0.7, jitter.width = 0.2),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.2) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    facet_grid(index_actual ~ index_counterfactual) +
    geom_text(data = df.text %>%
                drop_na(label),
              mapping = aes(x = x, y = y, label = as.character(label)),
              size = 6,
              position = position_dodge(width = .7),
              hjust = 0.5) +
    coord_cartesian(xlim = c(0.5, 2.5),
                    ylim = c(-5, 105),
                    expand = F,
                    clip = "off") +
    labs(y = ylabel) +
    theme(text = element_text(size = 20),
          legend.position = "none",
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.y = element_text(margin = margin(0, 10, 0, 0)),
          panel.spacing.y = unit(0.8, "cm"),
          plot.margin = unit(c(0.2, 0.2, .8, .2), "cm"),
          axis.title.x = element_blank(),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
  # print(p)
}

actual_counterfactual_threeballs_plot = function(x, ylabel) {
  p = ggplot(data = x, 
              mapping = aes(x = ball,
                            y = value,
                            group = model,
                            fill = colorindex)) +
    stat_summary(fun = "mean",
                 geom = "bar",
                 color = "black",
                 position = position_dodge(0.9),
                 width = 0.9) +
    stat_summary(fun.data = "mean_cl_boot",
                 geom = "linerange",
                 color = "black",
                 position = position_dodge(0.9)) +
    geom_point(data = x %>% 
                 mutate(value = ifelse(model == "rating", value, NA)),
               position = position_jitterdodge(jitter.width = 0.3,
                                               jitter.height = 0,
                                               dodge.width = 0.9),
               shape = 21,
               color = "black",
               size = 1,
               alpha = 0.15) +
    facet_wrap(~clip, ncol = 8) +
    geom_text(data = df.text, 
              mapping = aes(x = x,
                            y = y,
                            label = as.character(label)),
              size = 5, 
              position = position_dodge(width = .9), 
              hjust = 0.5, 
              fontface = "bold") +
    coord_cartesian(xlim = c(0.5, 2.5), 
                    ylim = c(-5, 120), 
                    expand = F,
                    clip = "off") +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    scale_fill_manual(values = c("red2", "green2", "black")) +
    scale_color_manual(values = c("red2", "green2", "black")) +
    labs(y = ylabel) +
    theme(text = element_text(size = 16),
          legend.position = "none",
          axis.title.x = element_blank(),
          panel.spacing.y = unit(0.3, "cm"),
          panel.grid = element_blank(),
          strip.background = element_blank(),
          strip.text = element_blank(),
          panel.border = element_rect(colour = "black", fill = NA))
}

2 Experiment 1: 2 balls

2.1 Clips

2.2 Read in data

# causal judgments
df.exp1.causal = read.csv("../../data/experiment1_causal.csv", 
                          header = T,
                          stringsAsFactors = F)

# counterfactual judgments
df.exp1.counterfactual = read.csv("../../data/experiment1_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# Information about each clip 
df.exp1.clipinfo = tibble(clip = 1:18,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 
                                             1, 1, 0, 1, 1, 1, 1, 1, 1),
                          outcome_counterfactual = c(0, 0, 0, 1, 1, 1, 0, 0, 1, 
                                                     0, 1, 1, 0, 0, 1, 0, 1, 1),
                          index_actual = rep(c("actual miss", 
                                               "actual close", 
                                               "actual hit"), each = 6),
                          index_counterfactual = rep(rep(c("counterfactual miss",
                                                           "counterfactual close", 
                                                           "counterfactual hit"),
                                                         each = 2),
                                                     3)) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp1.causal.long = df.exp1.causal$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

df.exp1.counterfactual.long = df.exp1.counterfactual$data %>%
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:18, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  arrange(participant, clip, rating) %>% 
  left_join(df.exp1.clipinfo,
            by = "clip")

2.2.1 Demographics

# counterfactual condition 
df.exp1.counterfactual %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
41 32 8.8 19 7.5 3.12
# causal condition 
df.exp1.causal %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time"),
           sep = ",") %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  select(gender:time) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
41 33 10.3 21 9.65 6.58

2.3 Read in model predictions

2.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "2ball")]

df.exp1.model = files %>%
  map_dfr(~ read.csv(str_c(path, .))) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp1.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp1.counterfactual.means = df.exp1.counterfactual.long %>%
  group_by(clip) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>%
  ungroup() %>%
  left_join(df.exp1.clipinfo, by = "clip")
`summarise()` ungrouping output (override with `.groups` argument)
# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp1.counterfactual.model = df.exp1.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, 
                   ~ sum((.x$prediction*100 - 
                            df.exp1.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

2.3.2 Causal condition

df.exp1.causal.model = df.exp1.counterfactual.long %>% 
  group_by(clip) %>% 
  summarize(empirical = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.counterfactual.model %>% 
              select(-c(noise, sse)),
            by = "clip") %>% 
  rename(simulation = prediction) %>% 
  mutate(empirical = ifelse(outcome_actual == 1, 100 - empirical, empirical),
         simulation = ifelse(outcome_actual == 1, 100 - simulation, simulation))
`summarise()` ungrouping output (override with `.groups` argument)

2.4 Counterfactual condition

2.4.1 Plots

set.seed(1)
df.plot = df.exp1.counterfactual.long %>% 
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.)/2)) %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", "actual close", "actual hit")),
         index_counterfactual = factor(index_counterfactual,
                                       levels = c("counterfactual miss",
                                                  "counterfactual close",
                                                  "counterfactual hit")),
         model = factor(model, levels = c("rating", "prediction"))) %>% 
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) #swap clips 11 and 12


df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(x = df.plot, ylabel = "counterfactual judgment")
print(p)

# ggsave("../../figures/plots/exp1_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

2.4.2 Stats

2.4.2.1 Model evaluation

df.exp1.counterfactual.means %>% 
  left_join(df.exp1.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise), 
            simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  print_table()
noise simulation_r simulation_rmse
1 0.97 10.95

2.5 Causal condition

2.5.1 Plots

set.seed(1)

model.name = c("empirical")
# model.name = c("simulation")

df.plot = df.exp1.causal.long %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp1.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, simulation, empirical),
               names_to = "model",
               values_to = "value") %>% 
  filter(model %in% c(model.name, "rating")) %>%
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(
    colorindex = ifelse(outcome_actual == 1, 2, 1),
    colorindex = ifelse(model != "rating", 3, colorindex),
    colorindex = as.factor(colorindex),
    index.actual = factor(index_actual, 
                          levels = c("actual miss", "actual close", "actual hit")),
    index.counterfactual = factor(index_counterfactual, 
                                  levels = c("counterfactual miss", 
                                             "counterfactual close", 
                                             "counterfactual hit"))) %>%
  mutate(clipindex = ifelse(clip == 11, 2, clipindex),
         clipindex = ifelse(clip == 12, 1, clipindex)) # swap clips 11 and 12

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(1:18, each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(df.plot, "causal judgment")
print(p)

# ggsave("../../figures/plots/exp1_causal_bars.pdf",
#        width = 8,
#        height = 6)

2.5.2 Stats

2.5.2.1 Model evaluation

df.exp1.causal.long %>% 
  mutate(rating = abs(rating)) %>% 
  group_by(clip,
           outcome_actual, 
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>% 
  summarize(rating = mean(rating)) %>% 
  ungroup() %>% 
  left_join(df.exp1.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>% 
  summarize(empirical_r = cor(rating, empirical),
            empirical_rmse = rmse(rating, empirical),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>% 
  print_table()
`summarise()` regrouping output by 'clip', 'outcome_actual', 'outcome_counterfactual', 'index_actual' (override with `.groups` argument)
empirical_r empirical_rmse simulation_r simulation_rmse
0.96 8.57 0.93 15.15

2.5.2.2 Bayesian regression

# random intercepts model 
fit.brm_exp1_causal = brm(formula = rating ~ 1 + index_actual + index_counterfactual + 
                            outcome_actual + (1 + index_actual + index_counterfactual + 
                                                outcome_actual | participant),
                          data = df.exp1.causal.long,
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp1_causal")

2.5.2.3 Specific hypotheses

df.exp1.causal.posterior = df.exp1.causal.long %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp1_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)
`summarise()` regrouping output by 'clip', 'index_actual', 'index_counterfactual' (override with `.groups` argument)
func_posterior_difference = function(data, clip1, clip2){
  data %>% 
    mutate(difference = .data[[clip1]] - .data[[clip2]]) %>% 
    pull(difference) %>% 
    mean_hdci() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])")) %>% 
    print()  
}

func_posterior_difference(data = df.exp1.causal.posterior,
                          clip1 = "8",
                          clip2 = "13")
            difference
1 (0.26 [-6.51, 7.31])

2.5.3 Tables

fit.brm_exp1_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  # print_table(format = "latex")
  print_table()
term estimate std.error lower 95% HDI upper 95% HDI
intercept -18.33 3.36 -23.81 -12.78
index_actualactualclose 4.31 3.40 -1.40 9.83
index_actualactualhit 4.06 4.84 -4.01 11.88
index_counterfactualcounterfactualclose -17.52 3.37 -22.96 -12.01
index_counterfactualcounterfactualhit -61.37 4.38 -68.54 -54.30
outcome_actual 92.36 5.11 84.04 100.96

3 Experiment 2: Brick & Teleport

3.1 Clips

3.2 Read in data

df.exp2.causal_first = read.csv("../../data/experiment2_causal_first.csv",
                                header = T,
                                stringsAsFactors = F)
df.exp2.counterfactual_first = read.csv("../../data/experiment2_counterfactual_first.csv",
                                        header = T,
                                        stringsAsFactors = F)

# Information about each clip

df.exp2.clipinfo = tibble(clip = 1:20,
                          outcome_actual = c(0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 
                                             0, 1, 1, 1, 1, 1, 1, 1, 0, 1),
                          outcome_counterfactual = c(0, 0, 1, 0, 1, 1, 0, 0, 1, 0, 
                                                     1, 1, 0, 0, 1, 0, 1, 1, 1, 1),
                          index_actual = c(rep(c("actual miss", 
                                                 "actual close",
                                                 "actual hit"),
                                               each = 6),
                                           rep("practice", 2)),
                          index_counterfactual = c(rep(rep(c("counterfactual miss",
                                                             "counterfactual close",
                                                             "counterfactual hit"),
                                                           each = 2), 
                                                       3),
                                                   rep("practice", 2))) %>% 
  mutate(index_actual = factor(index_actual, levels = c("actual miss", 
                                                        "actual close", 
                                                        "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss",
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

# Structure data 
df.exp2.causal_first.long = df.exp2.causal_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("causal", "counterfactual"), each = 20), max(participant)),
         condition = "causal_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

df.exp2.counterfactual_first.long = df.exp2.counterfactual_first$data %>% 
  enframe() %>% 
  mutate(participant = 1:n()) %>% 
  separate_rows(value, sep = ",") %>% 
  mutate(name = rep(c("clip", "rating"), n()/2),
         order = rep(rep(1:40, each = 2), max(participant)),
         value = as.numeric(value)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(question = rep(rep(c("counterfactual", "causal"), each = 20), max(participant)),
         condition = "counterfactual_first") %>% 
  left_join(df.exp2.clipinfo,
            by = "clip") %>%
  arrange(participant, question, clip) %>%
  select(participant, question, clip, rating, everything())

# combine data from both conditions
df.exp2.combined.long = df.exp2.causal_first.long %>%
  bind_rows(df.exp2.counterfactual_first.long %>% 
              mutate(participant = participant + 
                       max(df.exp2.causal_first.long$participant)))

3.2.1 Demographics

# counterfactual condition 
df.exp2.counterfactual_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time) %>% 
  bind_rows(df.exp2.causal_first %>% 
  separate(extra,
           into = c("gender", "age", "difficulty", "smoothness", "time", "condition"),
           sep = ",") %>% 
  select(gender:time)) %>% 
  mutate(gender = factor(gender, levels = c(1, 2), labels = c("female", "male"))) %>% 
  mutate_if(is.character, ~ as.numeric(.)) %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
82 35 12.2 45 21.25 5.11

3.3 Read in model predictions

3.3.1 Counterfactual condition

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "teleport")]

df.exp2.model = files %>%
  map_dfr(~ read.csv(str_c(path, .))) %>%
  set_names(c("clip", "noise", "prediction")) %>%
  left_join(df.exp2.clipinfo, by = "clip") %>%
  mutate(prediction = ifelse(outcome_actual == 1, 1 - prediction, prediction))

# calculate mean counterfactual judgments 
df.exp2.counterfactual.means = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  # mutate(rating = abs(rating)) %>% 
  group_by(clip) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.clipinfo, by = "clip")
`summarise()` ungrouping output (override with `.groups` argument)
# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp2.counterfactual.model = df.exp2.model %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data, 
                   ~ sum((.x$prediction*100 - df.exp2.counterfactual.means$rating) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

3.3.2 Causal condition

# Model predictions based on counterfactual judgments
df.exp2.causal.model = df.exp2.combined.long %>%
  filter(question == "counterfactual",
    clip <= 18) %>%
  group_by(clip, condition) %>%
  summarize(rating = mean(rating)) %>%
  pivot_wider(names_from = condition,
              values_from = rating) %>% 
  left_join(df.exp2.combined.long %>%
              filter(question == "counterfactual", clip <= 18) %>%
              group_by(clip) %>%
              summarize(combined = mean(rating)),
            by = "clip") %>% 
  left_join(df.exp2.counterfactual.model,
            by = "clip") %>% 
  select(-c(sse, noise)) %>% 
  rename(simulation = prediction) %>% 
  mutate_at(.vars = vars(simulation, causal_first, counterfactual_first, combined),
            .funs = list(~ ifelse(outcome_actual == 1, 100 - ., .)))
`summarise()` regrouping output by 'clip' (override with `.groups` argument)
`summarise()` ungrouping output (override with `.groups` argument)

3.4 Counterfactual condition

3.4.1 Plots

set.seed(1)

df.plot = df.exp2.combined.long %>%
  filter(question == "counterfactual",
         clip <= 18) %>%
  mutate(clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.counterfactual.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = ifelse(outcome_counterfactual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual,
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")),
         model = factor(model,
                             levels = c("rating", "prediction")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

p = actual_counterfactual_plot(df.plot, "counterfactual judgment")
print(p)

# ggsave("../../figures/plots/exp2_counterfactual_bars.pdf",
#        width = 8,
#        height = 6)

3.4.2 Stats

3.4.2.1 Model evaluation

df.exp2.counterfactual.means %>% 
  left_join(df.exp2.counterfactual.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual")) %>% 
  summarize(noise = unique(noise),
            simulation_r = cor(rating, prediction),
            simulation_rmse = rmse(rating, prediction)) %>% 
  print_table()
noise simulation_r simulation_rmse
0.9 0.96 13.46

3.4.2.2 Bayesian regressions

3.4.2.2.1 Check for potential order effects
# random intercepts & slopes model 
fit.brm_exp2_counterfactual = brm(formula = rating ~ 1 + condition * 
                                    (index_actual + index_counterfactual) + 
                                    (1 + index_actual + 
                                       index_counterfactual | participant),
                                  data = df.exp2.combined.long %>% 
                                    na.omit() %>% 
                                    filter(question == "counterfactual"),
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp2_counterfactual")

3.4.3 Tables

fit.brm_exp2_counterfactual %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  # print_table(format = "latex")
  print_table()
term estimate std.error lower 95% HDI upper 95% HDI
intercept 13.38 2.58 9.12 17.55
conditioncounterfactual_first -1.11 3.65 -6.99 4.82
index_actualactualclose -1.48 2.78 -6.02 3.11
index_actualactualhit -4.42 2.40 -8.32 -0.46
index_counterfactualcounterfactualclose 37.94 3.33 32.39 43.35
index_counterfactualcounterfactualhit 73.80 4.17 67.00 80.57
conditioncounterfactual_first:index_actualactualclose -0.55 3.84 -7.11 5.73
conditioncounterfactual_first:index_actualactualhit -0.85 3.41 -6.52 4.85
conditioncounterfactual_first:index_counterfactualcounterfactualclose 1.23 4.79 -6.62 9.18
conditioncounterfactual_first:index_counterfactualcounterfactualhit -0.02 5.92 -9.86 9.55

3.5 Causal condition

3.5.1 Plots

set.seed(1)

func_exp2_causal_plot = function(condition.name, model.name){

df.plot = df.exp2.combined.long %>%
  filter(question == "causal", clip <= 18) %>%
  filter(condition %in% condition.name) %>%
  mutate(rating = abs(rating),
         clipindex = rep(c(1, 2), nrow(.) / 2)) %>%
  left_join(df.exp2.causal.model,
            by = c("clip", 
                   "outcome_actual", 
                   "outcome_counterfactual", 
                   "index_actual", 
                   "index_counterfactual")) %>%
  pivot_longer(cols = c("rating", all_of(model.name)),
               names_to = "model",
               values_to = "value") %>% 
  mutate(model = factor(model, levels = c("rating", model.name))) %>%
  mutate(colorindex = ifelse(outcome_actual == 1, 2, 1),
         colorindex = ifelse(model != "rating", 3, colorindex),
         colorindex = as.factor(colorindex),
         index_actual = factor(index_actual, 
                               levels = c("actual miss", 
                                          "actual close", 
                                          "actual hit")),
         index_counterfactual = factor(index_counterfactual, 
                                       levels = c("counterfactual miss", 
                                                  "counterfactual close", 
                                                  "counterfactual hit")))

df.text = df.plot %>%
  expand(index_actual, index_counterfactual, clipindex, model) %>%
  mutate(label = rep(c("1 (A)", "2 (B)", "3", "4", "5 (A)", "6 (B)",
                       "7", "8 (C)", "9", "10", "11", "12 (C)",
                       "13 (D)", "14 (E)", "15", "16", "17 (D)", "18 (E)"), each = 2),
         label = ifelse(model != "rating", NA, label),
         y = -15,
         x = clipindex,
         colorindex = NA)

actual_counterfactual_plot(df.plot, "causal judgment")
}

plot.list = list(func_exp2_causal_plot(condition.name = "counterfactual_first",
                                       model.name = "combined"),
                 func_exp2_causal_plot(condition.name = "causal_first",
                                       model.name = "combined"))

# model.name = "combined"
# condition.name = "counterfactual_first"
# condition.name = "causal_first"
# 
# func_exp2_causal_plot(condition.name = condition.name,
#                       model.name = model.name),
# 
# ggsave(str_c("../../figures/plots/exp2_", condition.name ,"_causal_bars.pdf"),
#        width = 8,
#        height = 6)

3.5.2 Stats

3.5.2.1 Model evaluation

df.data = df.exp2.combined.long %>%
  filter(question == "causal",
         clip <= 18) %>%
  mutate(rating = abs(rating)) %>%
  group_by(condition,
           clip,
           outcome_actual,
           outcome_counterfactual,
           index_actual,
           index_counterfactual) %>%
  summarize(rating = mean(rating)) %>%
  ungroup() %>%
  left_join(df.exp2.causal.model,
            by = c("clip",
                   "outcome_actual",
                   "outcome_counterfactual",
                   "index_actual",
                   "index_counterfactual"))
`summarise()` regrouping output by 'condition', 'clip', 'outcome_actual', 'outcome_counterfactual', 'index_actual' (override with `.groups` argument)
df.data %>% 
  group_by(condition) %>% 
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined)) %>% 
  print_table()
`summarise()` ungrouping output (override with `.groups` argument)
condition empirical_r empirical_rmse
causal_first 0.92 13.78
counterfactual_first 0.99 5.27
df.data %>%
  summarize(empirical_r = cor(rating, combined),
            empirical_rmse = rmse(rating, combined),
            simulation_r = cor(rating, simulation),
            simulation_rmse = rmse(rating, simulation)) %>%
  print_table()
empirical_r empirical_rmse simulation_r simulation_rmse
0.95 10.43 0.92 19.09

3.5.2.2 Bayesian regression

# random intercepts + slopes model 
fit.brm_exp2_causal = brm(formula = rating ~ 1 + condition * 
                            (index_actual + index_counterfactual + outcome_actual) +
                            (1 + index_actual + index_counterfactual + 
                               outcome_actual | participant),
                          data = df.exp2.combined.long %>% 
                            filter(question == "causal",
                                   clip <= 18),
                          cores = 2,
                          seed = 1,
                          file = "cache/brm_exp2_causal")

3.5.2.3 Specific hypotheses

df.exp2.causal.posterior = df.exp2.combined.long %>% 
  filter(question == "causal",
         clip <= 18) %>% 
  group_by(clip, index_actual, index_counterfactual, outcome_actual, condition) %>% 
  summarize(value = mean(rating)) %>% 
  ungroup() %>% 
  add_fitted_draws(newdata = .,
                   model = fit.brm_exp2_causal,
                   re_formula = NA) %>% 
  ungroup() %>% 
  select(clip, condition, .value, .draw) %>% 
  pivot_wider(names_from = clip,
              values_from = .value)
`summarise()` regrouping output by 'clip', 'index_actual', 'index_counterfactual', 'outcome_actual' (override with `.groups` argument)
# difference between clips 13 and 14 versus 8 
df.exp2.causal.posterior %>% 
    mutate(difference = (`13`  + `14`) / 2 - `8`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))
               difference
1 (-0.55 [-12.09, 10.97])
# difference between clips 5 and 6 versus 11 
df.exp2.causal.posterior %>% 
    mutate(difference = (`5`  + `6`) / 2 - `11`) %>% 
    pull(difference) %>% 
    mean_hdi() %>% 
    mutate_if(is.numeric, ~round(., 2)) %>% 
    summarize(difference = str_c("(", y, " [", ymin, ", ", ymax, "])"))
             difference
1 (-5.48 [-13.8, 2.73])

3.5.3 Tables

3.5.3.1 Bayesian regression

fit.brm_exp2_causal %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  # print_table(format = "latex")
  print_table()
term estimate std.error lower 95% HDI upper 95% HDI
intercept -24.07 3.50 -29.80 -18.27
conditioncounterfactual_first 11.47 4.91 3.29 19.40
index_actualactualclose 8.17 3.44 2.39 13.62
index_actualactualhit 13.42 4.90 5.42 21.30
index_counterfactualcounterfactualclose -30.20 3.72 -36.41 -24.36
index_counterfactualcounterfactualhit -50.29 4.60 -57.94 -42.84
outcome_actual 98.53 5.40 89.53 107.57
conditioncounterfactual_first:index_actualactualclose -5.39 4.94 -13.38 2.80
conditioncounterfactual_first:index_actualactualhit -16.97 7.02 -28.61 -5.75
conditioncounterfactual_first:index_counterfactualcounterfactualclose -8.13 5.17 -16.51 0.35
conditioncounterfactual_first:index_counterfactualcounterfactualhit -25.51 6.44 -36.53 -15.18
conditioncounterfactual_first:outcome_actual 10.72 7.60 -1.63 23.34

4 Experiment 3: 3 balls

4.1 Clips

4.2 Read in data

# clipinfo
df.exp3.clipinfo = tibble(clip = 1:32,
                          outcome_both = c(0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 
                                           0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1),
                          outcome_a = c(0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1),
                          outcome_b = c(0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 
                                        0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1),
                          outcome_none = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                                           1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1),
                          clipindex = rep(1:2, 16))

# COUNTERFACTUAL JUDGMENTS 
df.exp3.counterfactual = read.csv("../../data/experiment3_counterfactual.csv",
                                  header = T,
                                  stringsAsFactors = F)

# demographics
df.exp3.counterfactual.demographics = df.exp3.counterfactual$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", "condition"),
                    n()/6),
         participant = rep(1:(n()/6), each = 6)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(condition = factor(condition, levels = 20:21, labels = c("A", "B")),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, condition, gender, age, time, smoothness, difficulty)

# judgments
df.exp3.counterfactual.long = df.exp3.counterfactual$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "rating"), n()/2),
         participant = rep(1:(n()/64), each = 64),
         order = rep(rep(1:32, each = 2), n()/64)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  left_join(df.exp3.counterfactual.demographics %>% 
              select(participant, ball = condition),
            by = "participant") %>% 
  select(participant, ball, clip, everything()) %>% 
  arrange(participant, clip)

# CAUSAL JUDGMENTS
df.exp3.causal = read.csv("../../data/experiment3_causal.csv",
                          header = T,
                          stringsAsFactors = F)

# demographics 
df.exp3.causal.demographics = df.exp3.causal$extra %>% 
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("gender", "age", "difficulty", "smoothness", "time", 
                      "counterbalance", "replayType", "experiment"),
                    n()/8),
         participant = rep(1:(n()/8), each = 8)) %>% 
  pivot_wider(names_from = name,
              values_from = value) %>% 
  mutate(counterbalance = ifelse(participant %in% c(1, 3, 5, 6, 9), 1, counterbalance),
         gender = factor(gender, levels = 1:2, labels = c("female", "male"))) %>% 
  select(participant, gender, age, counterbalance, time, smoothness, difficulty)

# judgments
df.exp3.causal.long = df.exp3.causal$data %>%
  enframe() %>% 
  separate_rows(value) %>% 
  mutate(value = as.numeric(value),
         name = rep(c("clip", "x", "y", "z", "rating1", "rating2"), n()/6),
         participant = rep(1:(n()/(32*6)), each = (32*6)),
         order = rep(rep(1:32, each = 6), n()/(32*6))) %>% 
  filter(!name %in% c("x", "y", "z")) %>% 
pivot_wider(names_from = name,
            values_from = value) %>% 
  left_join(df.exp3.causal.demographics %>% 
              select(participant, counterbalance),
            by = "participant") %>% 
  mutate(A = ifelse(counterbalance == 1, rating1, rating2),
         B = ifelse(counterbalance == 1, rating2, rating1)) %>% 
  select(-rating1, -rating2) %>% 
  pivot_longer(cols = c(A, B),
               names_to = "ball",
               values_to = "rating") %>% 
  select(participant, clip, ball, rating, order) %>% 
  arrange(participant, clip, ball)

4.2.1 Demographics

# counterfactual condition 
df.exp3.counterfactual.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
80 33 10.1 34 18.08 4.63
# causal condition 
df.exp3.causal.demographics %>% 
  summarize(n = nrow(.),
            age_mean = mean(age) %>% round(0),
            age_sd = sd(age) %>% round(1),
            n_female = sum(gender == "female"),
            time_mean = mean(time) %>% round(2),
            time_sd = sd(time) %>% round(2)) %>% 
  print_table()
n age_mean age_sd n_female time_mean time_sd
41 34 10.5 21 21.19 4.96

4.3 Read in model predictions

# read model predictions 
path = "../python/results/"
files = list.files(path = path, pattern = "*.csv")
files = files[str_detect(files, "3ball")]

df.exp3.model = files %>%
  map_dfr(~ read_csv(str_c(path, .))) %>% 
  rename(clip = trial)

4.3.1 Counterfactual condition

# calculate mean counterfactual judgments 
df.exp3.counterfactual.means = df.exp3.counterfactual.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()
`summarise()` regrouping output by 'clip' (override with `.groups` argument)
# find noisy simulation model that best predicts the mean counterfactual judgments by 
# minimizing the sum of squared errors 
df.exp3.model.counterfactual = df.exp3.model %>% 
  group_by(noise) %>% 
  select(clip, contains("whether"), noise) %>% 
  pivot_longer(cols = c(A_whether, B_whether),
               names_to = "ball",
               values_to = "prediction") %>% 
  mutate(ball = str_remove(ball, "_whether")) %>% 
  arrange(noise, clip, ball) %>% 
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  mutate(prediction = ifelse(outcome_both == 1, 1 - prediction, prediction)) %>% 
  group_by(noise) %>% 
  nest() %>% 
  mutate(sse = map(data,
                   ~ sum((.x$prediction*100 -
                            df.exp3.counterfactual.means$rating_mean) ^ 2)),
         sse = unlist(sse)) %>% 
  ungroup() %>% 
  filter(sse == min(sse)) %>% 
  unnest(data) %>% 
  mutate(prediction = prediction * 100)

4.3.2 Causal condition

# calculate mean causal judgments 
df.exp3.causal.means = df.exp3.causal.long %>% 
  group_by(clip, ball) %>%
  summarize(rating_mean = mean(rating),
            rating_low = smean.cl.boot(rating)[2],
            rating_high = smean.cl.boot(rating)[3]) %>% 
  ungroup()
`summarise()` regrouping output by 'clip' (override with `.groups` argument)
df.exp3.model.causal = df.exp3.model %>% 
  # take into account difference-making 
  mutate_at(.vars = vars(A_how:A_robust),
            .funs = list(~ . * A_difference)) %>% 
  mutate_at(.vars = vars(B_how:B_robust),
            .funs = list(~ . * B_difference)) %>% 
  # choose model based on best fit with counterfactual data 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>% 
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(clip, ball:robust)

4.3.3 Heuristic model

l.features = fromJSON("data/features.json")

df.features.initial = l.features[["appearance"]] %>%
  cbind(l.features[["initialVelocity"]][["ballA"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballB"]]) %>%
  cbind(l.features[["initialVelocity"]][["ballE"]]) %>%
  setNames(c(str_c("ball", LETTERS[c(1, 2, 5)], "_appearance"),
             "ballA_velx",
             "ballA_vely",
             "ballB_velx",
             "ballB_vely",
             "ballE_velx",
             "ballE_vely")) %>%
  mutate(clip = 1:nrow(.)) %>%
  pivot_longer(cols = starts_with("ball")) %>%
  separate(name, into = c("ball", "property")) %>%
  pivot_wider(names_from = property, values_from = value) %>% 
  mutate(ball = str_remove_all(ball, pattern = "ball"))

# information about collisions
df.features.collisions = data.frame()
ballnames = c("ballA", "ballB", "ballE")

for (i in 1:32) {
  df.tmp = data.frame()
  for (j in 1:length(ballnames)) {
    ncollisions = length(l.features[["collisions"]][[ballnames[j]]][["object"]][[i]])
    if (ncollisions > 0) {
      tmp = data.frame(
        ball = ballnames[j],
        object = l.features[["collisions"]][[ballnames[j]]][["object"]][[i]],
        time = l.features[["collisions"]][[ballnames[j]]][["time"]][[i]],
        pre.velx = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["x"]],
        pre.vely = l.features[["collisions"]][[ballnames[j]]][["pre"]][[i]][["y"]],
        post.velx = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["x"]],
        post.vely = l.features[["collisions"]][[ballnames[j]]][["post"]][[i]][["y"]],
        ncollision = 1:ncollisions
      )
      df.tmp = df.tmp %>%
        rbind(tmp)
    }
  }
  df.features.collisions = df.features.collisions %>%
    rbind(df.tmp %>%
            mutate(clip = i) %>%
            select(clip, ball, everything()))
}

# find end of clip
tmp = df.features.collisions %>%
  group_by(clip) %>%
  filter(ball == "ballE",
         str_detect(object, "goal|Left|Right")) %>%
  group_by(clip) %>%
  mutate(endclip = max(time)) %>%
  select(clip, endclip) %>%
  ungroup() %>%
  distinct()

# remove events that happen after the end of the clip
df.features.collisions = df.features.collisions %>%
  left_join(tmp, by = "clip") %>%
  mutate(endclip = ifelse(is.na(endclip), 400, endclip)) %>%
  group_by(clip) %>%
  filter(time <= endclip)

# E initially at rest or moving?
tmp.movingE = df.features.initial %>%
  filter(ball == "E") %>%
  mutate(E.moving = ifelse(velx == 0 & vely == 0, 1, 0)) %>%
  select(clip, E.moving)

# contact with ball E
tmp.contactE = df.features.collisions %>%
  filter(object == "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  mutate(contactE = 1) %>%
  select(clip, ball, contactE) %>%
  group_by(clip, ball) %>%
  summarize(contactE = any(contactE %>% as.logical()) * 1) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate(contactE = ifelse(is.na(contactE), 0, contactE))

# number of collisions
tmp.ncollisions = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  group_by(clip, ball) %>%
  count() %>%
  rename(ncollision = n)

# exclusivity (no contact with ball E by the other ball)
tmp.exclusive = df.features.collisions %>%
  filter(str_detect(object, "ball"),
         ball != "ballE") %>%
  count(clip, ball, object) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, 0),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, 0)) %>%
  group_by(clip) %>%
  summarize(AE = any(AE %>% as.logical()) * 1,
            BE = any(BE %>% as.logical()) * 1) %>%
  mutate(A = ifelse(AE == 1 & BE == 0, 1, 0),
         B = ifelse(AE == 0 & BE == 1, 1, 0)) %>%
  select(clip, A, B) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "exclusive")

# impact
func_angle = function(x, y) {
  dot.prod = x %*% y
  norm.x = norm(x, type = "2")
  norm.y = norm(y, type = "2")
  theta = acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}

# impact on ball E 
tmp.impactE = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ball"),
         ball == "ballE") %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(E.speed.diff = speed.diff,
         E.direction.diff = direction.diff)

# total impact 
tmp.impactTotal = df.features.collisions %>%
  rowwise() %>%
  filter(str_detect(object, "ballA|ballB")) %>%
  mutate(pre.speed = abs(pre.velx) + abs(pre.vely),
         post.speed = abs(post.velx) + abs(post.vely),
         speed.diff = post.speed - pre.speed,
         direction.diff = func_angle(c(pre.velx, pre.vely), c(post.velx, post.vely)),
         direction.diff = ifelse(is.na(direction.diff),
                                 abs(atan2(post.vely - pre.vely, post.velx - pre.velx)),
                                 direction.diff)) %>%
  ungroup() %>%
  group_by(clip, object) %>%
  summarize(speed.diff = sum(speed.diff),
            direction.diff = sum(direction.diff)) %>%
  rename(ball = object) %>%
  mutate(ball = factor(ball, levels = c("ballA", "ballB"), labels = c("A", "B"))) %>%
  left_join(expand.grid(clip = 1:32, ball = c("A", "B")), .,
            by = c("clip", "ball")) %>%
  mutate_at(.vars = vars(speed.diff, direction.diff),
            .funs = ~ ifelse(is.na(.), 0, .)) %>%
  arrange(clip, ball) %>%
  rename(total.speed.diff = speed.diff,
         total.direction.diff = direction.diff)

# transfer of force 
tmp.transfer = df.features.collisions %>%
  filter(str_detect(ball, "ballA|ballB"),
         str_detect(object, "ball")) %>%
  mutate(AE = ifelse(ball == "ballA" & object == "ballE", 1, NA),
         BE = ifelse(ball == "ballB" & object == "ballE", 1, NA),
         AB = ifelse((ball == "ballA" & object == "ballB"), 1, NA)) %>%
  mutate_at(.vars = vars(AE, BE, AB), .funs = ~ . * time) %>%
  filter(!is.na(AE) | !is.na(BE) | !is.na(AB)) %>%
  arrange(clip, time) %>%
  group_by(clip) %>%
  summarize(A = any(!is.na(AE)),
            B = any(!is.na(BE)),
            A = ifelse(any(!is.na(AB)) &
                         any(!is.na(BE)) &
                         max(AB, na.rm = T) < min(BE, na.rm = T),
                       T,
                       A),
            B = ifelse(any(!is.na(AB)) &
                         any(!is.na(AE)) &
                         max(AB, na.rm = T) < min(AE, na.rm = T),
                       T,
                       B)) %>%
  pivot_longer(cols = -clip,
               names_to = "ball",
               values_to = "transfer") %>% 
  arrange(clip, ball)

# collect features
df.features = df.features.initial %>%
  filter(ball != "E") %>%
  mutate(ball = factor(ball)) %>%
  mutate(moving = ifelse(velx == 0 & vely == 0, 0, 1),
         speed = abs(velx) + abs(vely)) %>%
  left_join(tmp.contactE) %>%
  left_join(tmp.impactE) %>%
  left_join(tmp.impactTotal) %>%
  left_join(tmp.transfer %>% mutate(transfer = transfer * 1)) %>%
  left_join(tmp.movingE) %>%
  left_join(tmp.exclusive) %>%
  select(-c(appearance, velx, vely))

# remove temporary variables
rm(list = ls()[str_detect(ls(), "tmp")])

4.4 Counterfactual condition

4.4.1 Plots

4.4.1.1 Bar graph

set.seed(1)

df.plot = df.exp3.counterfactual.long %>%
  left_join(df.exp3.model.counterfactual %>%
              select(clip, ball, prediction) %>%
              mutate(ball = as.factor(ball)),
            by = c("clip", "ball")) %>%
  left_join(df.exp3.clipinfo,
            by = "clip") %>% 
  pivot_longer(cols = c(rating, prediction),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" & ball == "A" & outcome_b == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "rating" & ball == "B" & outcome_a == 1,
                             2,
                             colorindex),
         colorindex = ifelse(model == "prediction", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "prediction"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>%
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

p = actual_counterfactual_threeballs_plot(df.plot, "counterfactual judgment")
print(p)
# ggsave("../../figures/plots/exp3_counterfactual_bars.pdf",
#        width = 12,
#        height = 6)

4.4.2 Stats

4.4.2.1 Model evaluation

# best fitting model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.model.counterfactual,
            by = c("clip", "ball")) %>% 
  summarize(simulation_r = cor(rating_mean, prediction),
            simulation_rmse = rmse(rating_mean, prediction)) %>% 
  print_table()
simulation_r simulation_rmse
0.86 19.84
# deterministic model 
df.exp3.counterfactual.means %>% 
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_a, outcome_b) %>% 
              pivot_longer(cols = -clip,
                           names_to = "ball",
                           values_to = "outcome") %>% 
              mutate(ball = factor(ball,
                                   levels = c("outcome_a", "outcome_b"),
                                   labels = c("B", "A"))),
            by = c("clip", "ball")) %>% 
  mutate(outcome = outcome * 100) %>% 
  summarize(simulation_r = cor(rating_mean, outcome),
            simulation_rmse = rmse(rating_mean, outcome)) %>% 
  print_table()
simulation_r simulation_rmse
0.83 30.28

4.5 Causal condition

4.5.1 Stats

4.5.1.1 Bayesian mixed effects model

df.data = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# using whether-causation from the model 
fit.brm_exp3_causal_w = brm(formula = rating ~ 1 + whether + (1 + whether | participant),
                                  data = df.data,
                                  save_all_pars = T,
                                  cores = 2,
                                  seed = 1,
                                  file = "cache/brm_exp3_causal_w")

fit.brm_exp3_causal_wh = brm(formula = rating ~ 1 + whether + how +
                               (1 + whether + how | participant),
                                   data = df.data,
                                   save_all_pars = T,
                                   cores = 2,
                                   seed = 1,
                                   file = "cache/brm_exp3_causal_wh")

fit.brm_exp3_causal_whs = brm(formula = rating ~ 1 + whether + how + sufficient +
                                (1 + whether + how + sufficient | participant),
                                    data = df.data,
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_whs")

# fit.brm_exp3_causal_w %>% summary()
# fit.brm_exp3_causal_wh %>% summary()
# fit.brm_exp3_causal_whs %>% summary()

4.5.1.2 Model comparison

fit.brm_exp3_causal_w = add_criterion(fit.brm_exp3_causal_w,
                                            criterion = c("loo", "waic"),
                                            reloo = T)
No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_wh = add_criterion(fit.brm_exp3_causal_wh,
                                             criterion = c("loo", "waic"),
                                             reloo = T)
No problematic observations found. Returning the original 'loo' object.
fit.brm_exp3_causal_whs = add_criterion(fit.brm_exp3_causal_whs,
                                              criterion = c("loo", "waic"),
                                              reloo = T)
No problematic observations found. Returning the original 'loo' object.
loo_compare(fit.brm_exp3_causal_w,
            fit.brm_exp3_causal_wh,
            fit.brm_exp3_causal_whs)
                        elpd_diff se_diff
fit.brm_exp3_causal_whs    0.0       0.0 
fit.brm_exp3_causal_wh  -125.2      15.6 
fit.brm_exp3_causal_w   -426.2      31.1 

4.5.1.3 Heuristic model

  • z-scored predictors and invidiual responses
# Regression based on features
df.regression.features = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo, by = "clip") %>% 
  left_join(df.features %>% 
              mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
            by = c("clip", "ball")) %>%
  clean_names() %>% 
  mutate(e_moving = 1 - e_moving)

# restrict the weights such that all predictors are positive
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

fit.brm_exp3_causal_heuristic = brm(formula = rating ~ moving + 
                                      speed + 
                                      contact_e + 
                                      e_speed_diff + 
                                      e_direction_diff + 
                                      total_speed_diff + 
                                      total_direction_diff + 
                                      transfer + 
                                      e_moving + 
                                      exclusive + (1 | participant),
                                    prior = priors,
                                    data = df.regression.features %>% 
                                      select(-c(clip, ball, order, clipindex,
                                                contains("outcome"))),
                                    save_all_pars = T,
                                    cores = 2,
                                    seed = 1,
                                    file = "cache/brm_exp3_causal_heuristic")
fit.brm_exp3_causal_heuristic
 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: rating ~ moving + speed + contact_e + e_speed_diff + e_direction_diff + total_speed_diff + total_direction_diff + transfer + e_moving + exclusive + (1 | participant) 
   Data: df.regression.features %>% select(-c(clip, ball, o (Number of observations: 2624) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~participant (Number of levels: 41) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     8.53      1.22     6.44    11.18 1.00     1163     1936

Population-Level Effects: 
                     Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept               49.79      1.46    46.91    52.74 1.00     1129
moving                   0.22      0.22     0.01     0.81 1.00     3399
speed                    2.06      0.85     0.44     3.67 1.00     2319
contact_e                0.38      0.35     0.01     1.29 1.00     3201
e_speed_diff             0.12      0.12     0.00     0.46 1.00     4435
e_direction_diff         1.04      0.74     0.04     2.71 1.00     2202
total_speed_diff         2.21      0.95     0.41     4.17 1.00     2208
total_direction_diff     3.99      0.91     2.14     5.68 1.00     2718
transfer                15.59      0.79    14.06    17.14 1.00     3555
e_moving                 0.18      0.17     0.01     0.62 1.00     3296
exclusive                4.39      0.73     2.97     5.81 1.00     3526
                     Tail_ESS
Intercept                1854
moving                   1547
speed                    1266
contact_e                2020
e_speed_diff             2268
e_direction_diff         1685
total_speed_diff         1393
total_direction_diff     1742
transfer                 2992
e_moving                 1738
exclusive                2538

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    33.21      0.46    32.33    34.10 1.00     5313     2845

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

4.5.1.4 Predictions for mean causal judgments

func_fit_data = function(fit){
  fit %>% 
    fitted(newdata = df.exp3.model.causal %>% 
             left_join(df.features %>% 
                         mutate_at(.vars = vars(moving:exclusive), .funs = ~ scale(.)[,]),
                       by = c("clip", "ball")) %>%
             clean_names() %>% 
             mutate(e_moving = 1 - e_moving),
           re_formula = NA) %>% 
    as_tibble() %>% 
    pull(Estimate)
}

df.exp3.causal.regression = df.exp3.causal.means %>% 
  mutate(w = func_fit_data(fit.brm_exp3_causal_w),
         wh = func_fit_data(fit.brm_exp3_causal_wh),
         whs = func_fit_data(fit.brm_exp3_causal_whs),
         heuristic = func_fit_data(fit.brm_exp3_causal_heuristic))

4.5.1.5 Model evaluation

prediction = "whs"

df.exp3.causal.regression %>% 
  summarize(simulation_r = cor(rating_mean, .[[prediction]]),
            simulation_rmse = rmse(rating_mean, .[[prediction]])) %>% 
  print_table()
simulation_r simulation_rmse
0.87 13.07

4.5.1.6 Individual participant regressions

df.fit = df.exp3.causal.long %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball"))

# priors 
priors = c(set_prior("normal(0,10)", class = "b", lb = 0))

# initial model fits (used for compilation)
fit.brm_exp3_causal_individual_baseline = 
  brm(formula = rating ~ 1,
      data = df.fit %>% 
        filter(participant == 1),
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_baseline"))

fit.brm_exp3_causal_individual_w = 
  brm(formula = rating ~ 1 + whether,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_w"))

fit.brm_exp3_causal_individual_wh = 
  brm(formula = rating ~ 1 + whether + how,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_wh"))

fit.brm_exp3_causal_individual_whs = 
  brm(formula = rating ~ 1 + whether + how + sufficient,
      data = df.fit %>% 
        filter(participant == 1),
      prior = priors,
      save_all_pars = T,
      cores = 2,
      seed = 1,
      file = str_c("cache/brm_exp3_causal_individual_whs"))

# update model fits for different participants
df.exp3.causal.individual_fit = df.fit %>%
  group_by(participant) %>%
  nest() %>%
  ungroup() %>% 
  # fit model for each participant
  mutate(fit_baseline = map(.x = data,
                            .f = ~ update(fit.brm_exp3_causal_individual_baseline,
                                          newdata = .x)),
         fit_w = map(.x = data,
                     .f = ~ update(fit.brm_exp3_causal_individual_w,
                                   newdata = .x)),
         fit_wh = map(.x = data,
                      .f = ~ update(fit.brm_exp3_causal_individual_wh,
                                    newdata = .x)),
         fit_whs = map(.x = data,
                       .f = ~ update(fit.brm_exp3_causal_individual_whs,
                                     newdata = .x))) %>% 
  mutate(fit_baseline = map(.x = fit_baseline,
                            ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_w = map(.x = fit_w,
                     ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_wh = map(.x = fit_wh,
                      ~ add_criterion(.x, criterion = c("loo", "waic"))),
         fit_whs = map(.x = fit_whs,
                       ~ add_criterion(.x, criterion = c("loo", "waic"))),
         r_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ cor(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         r_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ cor(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         r_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ cor(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         rmse_w = map2_dbl(.x = data,
                        .y = fit_w,
                        .f = ~ rmse(.x$rating, .y %>% 
                                     fitted() %>% 
                                     .[, 1])),
         rmse_wh = map2_dbl(.x = data,
                         .y = fit_wh,
                         .f = ~ rmse(.x$rating, .y %>% 
                                      fitted() %>% 
                                      .[, 1])),
         rmse_whs = map2_dbl(.x = data,
                          .y = fit_whs,
                          .f = ~ rmse(.x$rating, .y %>% 
                                       fitted() %>% 
                                       .[, 1])),
         model_comparison = pmap(.l = list(baseline = fit_baseline,
                                           w = fit_w,
                                           wh = fit_wh,
                                           whs = fit_whs),
                                 .f = ~ loo_compare(..1, ..2, ..3, ..4)),
         best_model = map_chr(.x = model_comparison,
                              .f = ~ rownames(.) %>% .[1]),
         best_model = factor(best_model,
                             levels = c("..1", "..2", "..3", "..4"),
                             labels = c("baseline", "w", "wh", "whs")))

save(list = c("df.exp3.causal.individual_fit"),
     file = "data/exp3_causal_individual_fit.RData")

4.5.1.7 Load individual participant regressions

load("data/exp3_causal_individual_fit.RData")

# count how many participants are best fit by the different models
df.exp3.causal.individual_fit %>% 
  count(best_model) %>% 
  print_table()
best_model n
wh 2
whs 39

4.5.1.8 Clusters of individual participants’ responses

set.seed(1)

df.cluster_whs = fit.brm_exp3_causal_whs %>% 
  ranef() %>% 
  .$participant %>% 
  as_tibble() %>% 
  select(contains("Estimate")) %>% 
  set_names(tolower(str_replace_all(names(.), "Estimate.", ""))) %>% 
  mutate(participant = 1:n()) %>% 
  relocate(participant)

df.cluster_whs = df.cluster_whs %>% 
  mutate(cluster = df.cluster_whs %>%
           select(-participant) %>%
           kmeans(centers = 2) %>% 
           .$cluster)

4.5.2 Plots

4.5.2.1 Bar graph

set.seed(1)
model_index = "whs"

# model predictions 
model_prediction = list(fit.brm_exp3_causal_w,
                        fit.brm_exp3_causal_wh,
                        fit.brm_exp3_causal_whs) %>%
  map_dfr(~ fitted(., df.exp3.causal.means %>% 
                     left_join(df.exp3.model.causal,
                               by = c("clip", "ball")),
                   re_formula = NA) %>% 
            as_tibble()) %>% 
  mutate(ball = rep(c("A", "B"), n()/2),
         clip = rep(rep(1:32, each = 2), 3),
         model = rep(c("w", "wh", "whs"),
                     each = 64))

df.plot = df.exp3.causal.long %>%
  left_join(df.exp3.clipinfo %>% 
              select(clip, outcome_both),
            by = c("clip")) %>% 
  left_join(model_prediction %>% 
              filter(model == model_index) %>% 
              select(model = Estimate, clip, ball),
            by = c("clip", "ball")) %>% 
  pivot_longer(cols = c(rating, model),
               names_to = "model",
               values_to = "value") %>% 
  mutate(colorindex = 1,
         colorindex = ifelse(model == "rating" & outcome_both == 0, 1, colorindex),
         colorindex = ifelse(model == "rating" & outcome_both == 1, 2, colorindex),
         colorindex = ifelse(model == "model", 3, colorindex),
         colorindex = as.factor(colorindex),
         model = factor(model, levels = c("rating", "model"))) %>%
  arrange(participant, clip, ball)

df.text = df.plot %>%
  expand(ball, clip, model) %>% 
  mutate(label = ifelse(model == "rating" & ball == "A", clip, NA),
         y = 110,
         x = 1.5,
         colorindex = NA)

p = actual_counterfactual_threeballs_plot(df.plot,"causal responsibility")
print(p)

# ggsave("../../figures/plots/exp3_causal_bars.pdf",
#        width = 12,
#        height = 6)

4.5.2.2 Selection

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))
  

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(3, 7, 15, 16, 23)) %>%
  group_by(clip, ball) %>%
  summarize(mean = mean(rating),
            low = smean.cl.boot(rating)[2],
            high = smean.cl.boot(rating)[3]) %>%
  left_join(df.exp3.causal.regression %>% select(clip, ball, w, wh, whs),
            by = c("clip", "ball")) %>%
  ungroup() %>%
  pivot_longer(cols = c(mean, w, wh, whs),
               names_to = "index",
               values_to = "value") %>% 
  mutate_at(.vars = vars(low, high), .funs = ~ ifelse(index == "mean", ., NA)) %>%
  mutate(index = factor(index, levels = c("mean", "w", "wh", "whs")),
         clip = factor(clip, levels = c(7, 23, 3, 15, 16),
                       labels = c("causal chain", "double prevention", "joint causation",
                                  "overdetermination", "preemption")))
`summarise()` regrouping output by 'clip' (override with `.groups` argument)
df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         value = -10,
         index = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         index = NA,
         value = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = value,
                     group = index,
                     fill = index)) +
  geom_bar(stat = "identity", color = "black",
           position = position_dodge(0.9),
           width = 0.9) +
  geom_errorbar(mapping = aes(ymin = low, ymax = high),
                width = 0,
                size = 1,
                position = position_dodge(0.9)) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -25, ymax = -7) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -25, ymax = -7) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) + 
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_grid(~clip) +
  scale_fill_grey(start = 1,
                  end = 0,
                  labels = c("mean rating  ", expression(CSM[W] ~ " "),
                             expression(CSM[WH] ~ " "),
                             expression(CSM[WHS] ~ " "))) +
  labs(y = "causal responsibility", fill = "") +
  coord_cartesian(clip = "off") + 
  theme_bw() +
  theme(legend.text.align = 0,
        text = element_text(size = 20),
        panel.grid = element_blank(),
        legend.position = "bottom",
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        legend.text = element_text(size = 30),
        strip.text = element_text(size = 30),
        legend.spacing = unit(0.5, "cm"),
        legend.background = element_rect(fill = "transparent"),
        legend.margin = margin(t = 5, unit = "cm"),
        plot.margin = margin(t = 8, l = 0.2, r = 0.2, b = 0.1, unit = "cm"))
print(p)

# ggsave("../../figures/plots/exp3_selection_bars.pdf",
#        width = 20,
#        height = 10,
#        device = cairo_pdf)

4.5.2.3 Scatter plot

func_scatterplot = function(model){
  if(model == "heuristic"){
    xlabel = "Heuristic"
  }else{
    xlabel = bquote(CSM[.(toupper(model))])
  }
  
  l.models = list(w = fit.brm_exp3_causal_w, 
                  wh = fit.brm_exp3_causal_wh,
                  whs = fit.brm_exp3_causal_whs,
                  heuristic = fit.brm_exp3_causal_heuristic)
  
  df.data = df.exp3.causal.means %>% 
    left_join(df.exp3.model.causal, by = c("clip", "ball")) %>% 
    left_join(df.regression.features, by = c("clip", "ball"))
  
  df.plot = l.models[[model]] %>%
    fitted(newdata = df.data,
           re_formula = NA) %>% 
    as_tibble() %>% 
    clean_names() %>%
    bind_cols(df.data)
  
  p = ggplot(data = df.plot,
             mapping = aes(x = estimate, 
                           y = rating_mean)) +
    geom_abline(intercept = 0,
                slope = 1,
                linetype = 2) + 
    geom_smooth(aes(y = estimate,
                    ymin = q2_5,
                    ymax = q97_5),
                stat = "identity",
                color = "black") +
    geom_linerange(size = 0.75, 
                   mapping = aes(ymin = rating_low,
                                 ymax = rating_high),
                   color = "gray50") +
    geom_point(size = 2) + 
    scale_color_manual(values = c("black", "#e41a1c"),
                       guide = F) +
    coord_fixed(xlim = c(0, 100),
                ylim = c(0, 100)) +
    labs(x = xlabel,
         y = "responsibility judgment") +
    annotate(geom = "text",
             label = str_c(
               "r = ", cor(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2) %>% 
                 as.character() %>% 
                 str_sub(start = 2),
               "\nRMSE = ", rmse(df.plot$estimate, df.plot$rating_mean) %>% 
                 round(2)),
             x = 5,
             y = 95,
             size = 6,
             hjust = 0) +
    scale_x_continuous(breaks = seq(0, 100, 25)) +
    scale_y_continuous(breaks = seq(0, 100, 25)) +
    theme(text = element_text(size = 20))
  return(p)
}

# for creating and saving an individual scatter plot 
# model = "w"
# model = "wh"
# model = "whs"
# model = "heuristic"

# func_scatterplot(model)
# ggsave(str_c("../../figures/plots/exp3_", model, "_scatter.pdf"),
#        width = 5,
#        height = 5)

list(func_scatterplot(model = "w"),
     func_scatterplot(model = "wh"),
     func_scatterplot(model = "whs"),
     func_scatterplot(model = "heuristic")) %>%
  wrap_plots(ncol = 2)

4.5.2.4 Individual participant variance for a selection of clips (clustered)

# \u2713 = check mark 
# \u2717 = cross mark 

check = "\u2713"
cross = "\u2717"

x_labels = c(str_c("\nW ", cross, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", cross, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", check, "\nH ", check, "\nS ", cross),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", check, "\nS ", check),
             str_c("\nW ", cross, "\nH ", cross, "\nS ", cross))

df.plot = df.exp3.causal.long %>%
  filter(clip %in% c(7, 23, 3, 15, 16)) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 3, 15, 16),
    labels = c("causal chain", "double prevention", "joint causation",
      "overdetermination", "preemption")))

df.plot = df.plot %>% 
  left_join(df.cluster_whs %>% 
              group_by(cluster) %>% 
              mutate(group = factor(cluster, levels = 1:2, labels = str_c("n = ", group_size(.)))) %>% 
              ungroup() %>% 
              select(participant, cluster, group),
            by = "participant")


df.labels = df.plot %>% 
  distinct(clip, ball) %>% 
  arrange(clip, ball) %>% 
  mutate(labels = x_labels,
         rating = -10,
         group = NA,
         participant = NA)

func_load_image = function(clip){
  readPNG(str_c("../../figures/diagrams/exp3_clip", clip, ".png"))
}

df.clips = df.plot %>% 
  distinct(clip) %>% 
  arrange(clip) %>% 
  mutate(number = c(7, 23, 3, 15, 16),
         grob = map(.x = number, .f = ~ func_load_image(clip = .x)),
         group = NA,
         participant = NA,
         ball = NA,
         label = str_c("clip ", number))

ball_a = readPNG("../../figures/diagrams/ballA.png")
ball_a = rasterGrob(ball_a, interpolate = TRUE)

ball_b = readPNG("../../figures/diagrams/ballB.png")
ball_b = rasterGrob(ball_b, interpolate = TRUE)

p = ggplot(data = df.plot, 
       mapping = aes(x = ball,
                     y = rating,
                     group = participant,
                     shape = group)) +
  geom_line(mapping = aes(linetype = "individual",
                          color = group),
            size = 1,
            alpha = 0.3) +
  geom_point(mapping = aes(color = group),
             size = 1,
             alpha = 0.3) +
  stat_summary(mapping = aes(group = group,
                             color = group,
                             linetype = "mean"),
               fun = "mean",
               geom = "line",
               size = 1.5) +
  stat_summary(data = df.plot %>% filter(cluster == 1),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  stat_summary(data = df.plot %>% filter(cluster == 2),
               mapping = aes(group = group,
                             color = group),
               fun.data = "mean_cl_boot",
               geom = "pointrange",
               size = 1,
               shape = 19) +
  annotation_custom(grob = ball_a, xmin = 0.5, xmax = 1.5, ymin = -30, ymax = -10) + 
  annotation_custom(grob = ball_b, xmin = 1.5, xmax = 2.5, ymin = -30, ymax = -10) +
  geom_text(data = df.labels,
            mapping = aes(label = labels,
                          y = -Inf),
            vjust = 1.2,
            family = "Arial Unicode MS",
            size = 8) + 
  geom_custom(data = df.clips,
              mapping = aes(data = grob, x = 1.5, y = Inf),
              grob_fun = function(x) rasterGrob(x,
                                                interpolate = T,
                                                vjust = -0.25)) +
  geom_text(data = df.clips,
            mapping = aes(label = label,
                          y = Inf,
                          x = -Inf),
            size = 9,
            hjust = -0.2,
            vjust = -4) + 
  facet_wrap(~ clip,
             ncol = 8) +
  coord_cartesian(xlim = c(0.9, 2.1),
                  ylim = c(-5, 105),
                  expand = F,
                  clip = "off") +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     labels = seq(0, 100, 25)) +
  scale_color_brewer(palette = "Set1",
                     guide = "none") +
  labs(y = "causal responsibility rating",
       linetype = "legend",
       shape = "") +
  theme_bw() +
  guides(linetype = guide_legend(override.aes = list(alpha = c(0.3, 1)),
                                 keywidth = unit(1.2, "cm")),
         shape = guide_legend(override.aes = list(color = c("#377eb8", "#e41a1c"),
                                                  shape = c(19, 19),
                                                  alpha = c(1, 1),
                                                  size = c(5, 5)))) +
  theme(panel.grid = element_blank(),
        text = element_text(size = 20),
        legend.position = c(0.505, 0.23),
        axis.title.x = element_blank(),
        legend.text = element_text(size = 20),
        legend.title = element_text(size = 24),
        legend.background = element_rect(fill = "transparent"),
        strip.text = element_text(size = 30),
        axis.text.x = element_blank(),
        legend.key = element_rect(fill = "transparent"),
        legend.box = "horizontal",
        legend.spacing.x = unit(0.1, "cm"),
        panel.spacing.x = unit(1, "cm"),
        plot.margin = margin(b = 5, l = 0.2, r = 0.2, t = 7.5, unit = "cm"))

print(p)

# ggsave(str_c("../../figures/plots/exp3_individual_variance_selection_lines_clustered.pdf"),
#        plot = p,
#        width = 20,
#        height = 8.5,
#        device = cairo_pdf)

4.5.3 Tables

4.5.3.1 Relationship between predictors

df.exp3.model %>% 
  filter(noise == unique(df.exp3.model.counterfactual$noise)) %>%
  pivot_longer(cols = A_difference:B_robust,
               names_to = c("ball", "aspect"),
               names_sep = "_",
               values_to = "value") %>% 
  pivot_wider(names_from = aspect,
              values_from = value) %>% 
  select(difference, whether, how, sufficient, robust) %>% 
  correlate() %>% 
  shave() %>% 
  fashion() %>% 
  # print_table(format = "latex")
  print_table()
rowname difference whether how sufficient robust
difference
whether .50
how .79 .27
sufficient .21 .10 .36
robust .43 .93 .24 .20

4.5.3.2 Population level predictors in the mixed effects models

fit.brm_exp3_causal_w %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(model = "CSM_w") %>% 
  bind_rows(fit.brm_exp3_causal_wh %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_wh")) %>% 
  bind_rows(fit.brm_exp3_causal_whs %>% 
              tidy() %>% 
              filter(str_detect(term, "b_")) %>% 
              mutate(model = "CSM_whs")) %>% 
  mutate(term = tolower(term),
         term = str_remove_all(term, "b_")) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  select(model, everything()) %>% 
  # print_table(format = "latex")
  print_table()
model term estimate std.error lower 95% HDI upper 95% HDI
CSM_w intercept 31.68 1.61 29.02 34.37
CSM_w whether 39.30 2.16 35.76 42.82
CSM_wh intercept 10.29 1.56 7.80 12.83
CSM_wh whether 28.40 2.68 23.96 32.70
CSM_wh how 36.07 2.60 31.80 40.27
CSM_whs intercept 10.43 1.59 7.82 13.01
CSM_whs whether 27.93 2.60 23.67 32.18
CSM_whs how 26.47 2.91 21.76 31.35
CSM_whs sufficient 30.90 3.06 25.77 35.96

4.5.3.3 Individual participants regression results

df.exp3.causal.individual_fit %>% 
  select_if(~ (is.numeric(.) | is.factor(.))) %>% 
  mutate_at(.vars = vars(contains("_w")), .funs = ~ round(., 2)) %>% 
  select(participant, everything(), best_model) %>% 
  # print_table(format = "latex")
  print_table()
participant r_w r_wh r_whs rmse_w rmse_wh rmse_whs best_model
1 0.29 0.38 0.48 30.00 29.03 27.74 whs
2 0.20 0.77 0.77 39.32 27.94 27.55 whs
3 0.37 0.78 0.79 38.90 28.00 26.95 whs
4 0.41 0.55 0.63 43.62 40.60 38.31 whs
5 0.45 0.53 0.54 39.23 37.32 36.81 whs
6 0.18 0.54 0.54 40.84 36.32 36.07 whs
7 0.49 0.61 0.64 32.96 29.93 29.15 whs
8 0.39 0.63 0.68 38.37 33.26 31.64 whs
9 0.40 0.72 0.76 35.86 28.34 26.39 whs
10 0.28 0.52 0.53 29.13 26.28 25.85 whs
11 0.51 0.56 0.60 33.66 32.19 31.17 whs
12 0.56 0.57 0.70 32.78 32.08 28.84 whs
13 0.56 0.71 0.74 29.83 25.40 24.21 whs
14 0.43 0.47 0.50 33.72 32.82 32.23 whs
15 0.19 0.35 0.37 39.77 38.29 37.94 whs
16 0.62 0.67 0.76 40.49 37.98 34.39 whs
17 0.42 0.66 0.71 28.57 23.82 22.42 whs
18 0.31 0.35 0.39 37.78 37.09 36.57 whs
19 0.37 0.60 0.62 33.41 29.23 28.62 whs
20 0.57 0.61 0.68 36.47 35.00 32.71 whs
21 0.43 0.50 0.57 30.95 29.61 28.32 whs
22 0.51 0.67 0.71 33.72 29.36 28.00 whs
23 0.33 0.45 0.52 38.76 37.00 35.57 whs
24 0.46 0.63 0.69 34.01 29.96 28.22 whs
25 0.61 0.66 0.70 30.60 28.91 27.55 whs
26 0.57 0.71 0.76 40.43 35.03 32.80 whs
27 0.46 0.52 0.66 43.12 41.49 38.12 whs
28 0.30 0.45 0.46 29.33 27.62 27.38 whs
29 0.48 0.60 0.65 38.93 36.03 34.48 whs
30 0.18 0.70 0.70 33.06 24.95 24.71 whs
31 0.48 0.59 0.69 37.54 34.72 31.76 whs
32 0.37 0.61 0.70 40.62 35.94 32.92 whs
33 0.50 0.56 0.64 38.64 36.79 34.85 whs
34 0.24 0.51 0.52 44.50 40.86 40.19 whs
35 0.15 0.32 0.30 34.56 33.36 33.37 wh
36 0.19 0.84 0.82 45.38 29.58 29.70 wh
37 0.55 0.63 0.67 36.04 33.36 32.11 whs
38 0.52 0.61 0.67 36.03 33.39 31.53 whs
39 0.47 0.77 0.77 29.61 21.75 21.58 whs
40 0.46 0.78 0.79 32.41 23.60 22.69 whs
41 0.32 0.50 0.58 32.90 30.38 28.69 whs

4.5.3.4 CSMwhs predictions for selection of cases

df.exp3.model.causal %>%
  left_join(df.exp3.causal.regression,
            by = c("clip", "ball")) %>% 
  filter(clip %in% c(7, 23, 16, 3, 15)) %>%
  mutate_at(.vars = vars(difference, robust, sufficient, whether, heuristic),
            .funs = ~ round(., 2)) %>%
  select(clip, ball, difference, whether, how, sufficient, robust) %>%
  mutate(clip = factor(clip, levels = c(7, 23, 16, 3, 15))) %>%
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ ifelse(. < 0.5,
                             str_c("xmark (", ., ")"),
                             str_c("cmark (", ., ")"))) %>%
  arrange(clip) %>% 
  # print_table(format = "latex")
  print_table()
clip ball difference whether how sufficient robust
7 A cmark (1) xmark (0.34) cmark (1) xmark (0) xmark (0.25)
7 B cmark (1) cmark (1) cmark (1) cmark (0.67) cmark (0.6)
23 A xmark (0.05) xmark (0) xmark (0) xmark (0) xmark (0)
23 B cmark (0.91) cmark (0.79) xmark (0) xmark (0) cmark (0.72)
16 A cmark (1) xmark (0.23) cmark (1) cmark (1) xmark (0.35)
16 B xmark (0) xmark (0) xmark (0) xmark (0) xmark (0)
3 A cmark (1) cmark (0.88) cmark (1) xmark (0.12) cmark (0.76)
3 B cmark (1) cmark (0.89) cmark (1) xmark (0.11) cmark (0.75)
15 A cmark (1) xmark (0.01) cmark (1) cmark (0.99) xmark (0.1)
15 B cmark (1) xmark (0.01) cmark (1) cmark (1) xmark (0.1)

4.5.3.5 CSMwhs predictions for all cases

df.exp3.causal.regression %>% 
  left_join(df.exp3.model.causal,
            by = c("clip", "ball")) %>% 
  left_join(df.exp3.clipinfo %>% 
              select(-clipindex),
            by = c("clip")) %>% 
  mutate_at(.vars = vars(difference, whether, how, sufficient, robust),
            .funs = ~ . * 100) %>% 
  mutate_at(.vars = vars(-ball), .funs = ~ round(., 0)) %>%
  select(clip, ball, contains("outcome"), difference, whether, how, sufficient, robust,
         w, wh, whs, heuristic, rating = rating_mean) %>% 
  # print_table(format = "latex", digits = 0)
  print_table(digits = 0)
clip ball outcome_both outcome_a outcome_b outcome_none difference whether how sufficient robust w wh whs heuristic rating
1 A 0 0 0 0 100 40 100 23 36 47 58 55 57 42
1 B 0 0 0 0 100 15 100 16 9 38 51 46 54 37
2 A 0 0 0 0 57 12 0 0 10 36 14 14 25 21
2 B 0 0 0 0 18 0 0 0 0 32 10 10 24 19
3 A 1 0 0 0 100 88 100 12 76 66 71 65 72 76
3 B 1 0 0 0 100 89 100 11 75 67 72 65 72 75
4 A 1 0 0 0 100 78 100 4 78 62 69 60 58 63
4 B 1 0 0 0 100 95 100 15 57 69 73 68 54 78
5 A 0 0 1 0 100 90 100 0 47 67 72 62 47 71
5 B 0 0 1 0 100 0 100 0 0 32 46 37 68 22
6 A 0 0 1 0 100 59 100 16 35 55 63 59 53 73
6 B 0 0 1 0 100 18 100 6 14 39 52 44 53 22
7 A 1 0 1 0 100 34 100 0 25 45 56 46 70 59
7 B 1 0 1 0 100 100 100 67 60 71 75 86 64 79
8 A 1 0 1 0 0 0 0 0 0 32 10 10 25 7
8 B 1 0 1 0 100 100 100 100 100 71 75 96 84 92
9 A 0 1 0 0 0 0 0 0 0 32 10 10 14 8
9 B 0 1 0 0 100 100 100 0 100 71 75 65 78 90
10 A 0 1 0 0 77 18 0 0 22 39 15 16 15 23
10 B 0 1 0 0 98 79 0 0 63 63 33 33 21 55
11 A 1 1 0 0 100 70 100 77 68 59 66 80 71 93
11 B 1 1 0 0 0 0 0 0 0 32 10 10 16 4
12 A 1 1 0 0 100 82 100 74 83 64 70 83 57 77
12 B 1 1 0 0 100 0 100 12 24 32 46 41 53 37
13 A 0 1 1 0 67 34 0 0 35 45 20 20 14 8
13 B 0 1 1 0 70 35 0 0 35 46 20 20 21 64
14 A 0 1 1 0 97 91 0 0 59 67 36 36 21 22
14 B 0 1 1 0 91 77 0 0 51 62 32 32 20 18
15 A 1 1 1 0 100 1 100 99 10 32 47 68 71 76
15 B 1 1 1 0 100 1 100 100 10 32 47 68 71 76
16 A 1 1 1 0 100 23 100 100 35 41 53 74 80 92
16 B 1 1 1 0 0 0 0 0 0 32 10 10 14 4
17 A 0 0 0 1 100 19 100 37 18 39 52 54 66 69
17 B 0 0 0 1 100 0 100 36 17 32 46 48 65 46
18 A 0 0 0 1 100 11 100 40 17 36 49 52 55 63
18 B 0 0 0 1 100 7 100 37 9 34 48 50 56 66
19 A 1 0 0 1 100 74 100 7 65 61 67 60 55 53
19 B 1 0 0 1 100 72 100 7 65 60 67 59 55 49
20 A 1 0 0 1 100 92 100 8 72 68 72 65 57 41
20 B 1 0 0 1 100 88 100 4 53 66 71 63 56 71
21 A 0 0 1 1 100 47 100 40 45 50 60 63 58 80
21 B 0 0 1 1 100 9 100 21 10 35 49 46 59 18
22 A 0 0 1 1 100 100 100 89 83 71 75 92 48 60
22 B 0 0 1 1 100 8 100 0 15 35 49 39 53 42
23 A 1 0 1 1 5 0 0 0 0 32 10 10 15 3
23 B 1 0 1 1 91 79 0 0 72 63 33 33 22 39
24 A 1 0 1 1 100 66 100 4 63 58 65 57 57 44
24 B 1 0 1 1 100 94 100 22 79 69 73 70 54 73
25 A 0 1 0 1 100 25 100 21 26 42 53 50 69 43
25 B 0 1 0 1 100 74 100 54 65 61 67 74 56 73
26 A 0 1 0 1 100 6 100 3 9 34 48 40 60 39
26 B 0 1 0 1 100 87 100 35 54 66 71 72 46 69
27 A 1 1 0 1 100 97 100 52 97 70 74 80 67 80
27 B 1 1 0 1 0 0 0 0 0 32 10 10 17 6
28 A 1 1 0 1 100 90 100 22 80 67 72 69 79 89
28 B 1 1 0 1 0 0 0 0 0 32 10 10 12 5
29 A 0 1 1 1 100 58 100 24 44 55 63 61 66 47
29 B 0 1 1 1 100 63 100 24 38 57 64 62 54 67
30 A 0 1 1 1 100 57 100 29 49 54 63 62 63 58
30 B 0 1 1 1 100 46 100 24 39 50 60 57 63 56
31 A 1 1 1 1 100 2 100 4 3 32 47 39 63 44
31 B 1 1 1 1 100 4 100 4 4 33 47 39 53 46
32 A 1 1 1 1 0 0 0 0 0 32 10 10 16 5
32 B 1 1 1 1 100 75 100 66 73 61 68 78 65 71

4.5.3.6 Heuristic model

fit.brm_exp3_causal_heuristic %>% 
  tidy() %>% 
  filter(str_detect(term, "b_")) %>% 
  mutate(term = str_remove(term, "b_"),
         term = tolower(term)) %>% 
  mutate_if(is.numeric, ~ round(., 2)) %>% 
  rename(`lower 95% HDI` = lower,
         `upper 95% HDI` = upper) %>% 
  # print_table(format = "latex")
  print_table()
term estimate std.error lower 95% HDI upper 95% HDI
intercept 49.79 1.46 47.39 52.23
moving 0.22 0.22 0.01 0.66
speed 2.06 0.85 0.64 3.43
contact_e 0.38 0.35 0.02 1.07
e_speed_diff 0.12 0.12 0.01 0.37
e_direction_diff 1.04 0.74 0.08 2.43
total_speed_diff 2.21 0.95 0.67 3.80
total_direction_diff 3.99 0.91 2.45 5.45
transfer 15.59 0.79 14.30 16.90
e_moving 0.18 0.17 0.01 0.51
exclusive 4.39 0.73 3.16 5.56

5 Appendix

5.0.0.1 Experiment 3

5.0.0.1.1 Does the question framing matter? (Responsibility vs. causation)

We have run a pilot eye-tracking experiment with 15 participants who saw the same clips as the ones in Experiment 3. Instead of asking participants to judge “To what extent were A and B responsible for E (not) going through the gate?”, we asked them to judge to what extent they agreed with the following statement: “Ball A/B caused ball E to go through the gate.” or “Ball A/B prevented ball E from going through the gate.” depending on the outcome.

Participants agreement judgments with the causal statements in this pilot eye-tracking experiment were strikingly similar to participants’ responsibility judgments in the online experiment.

df.exp3.eye_tracking_pilot = read_csv("../../data/experiment3_eye_tracking_pilot.csv") %>% 
  rename(clip = trial) %>%
  mutate(ball = str_to_upper(ball)) %>% 
  group_by(clip, ball, outcome) %>% 
  summarize(causality_mean = mean(causation),
            causality_low = smean.cl.boot(causation)[2],
            causality_high = smean.cl.boot(causation)[3]) %>% 
  ungroup()

df.plot = df.exp3.eye_tracking_pilot %>% 
  # select(clip, ball, outcome,) %>% 
  left_join(df.exp3.causal.means %>% 
              select(clip,
                     ball,
                     responsibility_mean = rating_mean,
                     responsibility_low = rating_low,
                     responsibility_high = rating_high),
            by = c("clip", "ball"))

ggplot(data = df.plot,
       mapping = aes(x = causality_mean,
                     y = responsibility_mean)) + 
  geom_abline(intercept = 0, 
              slope = 1,
              linetype = 2) + 
  geom_linerange(mapping = aes(ymin = responsibility_low,
                               ymax = responsibility_high),
                 alpha = 0.1) +
  geom_linerange(mapping = aes(xmin = causality_low,
                               xmax = causality_high),
                 alpha = 0.1) +
  geom_point(size = 2) + 
  geom_smooth(method = "lm",
              color = "black") +
  annotate(geom = "text",
           x = c(0, 0),
           y = c(100, 90),
           label = c(str_c("r = ",
                           cor(df.plot$causality_mean,df.plot$responsibility_mean) %>% 
                             round(2)),
                     str_c("RMSE = ", rmse(df.plot$causality_mean,df.plot$responsibility_mean) %>% 
                             round(2))),
           hjust = 0, 
           size = 8) + 
  scale_x_continuous(breaks = seq(0, 100, 25),
                     limits = c(0, 100))  +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     limits = c(0, 100))  +
  labs(x = "causal judgment",
       y = "responsibility judgment") +
  theme(text = element_text(size = 24))

# ggsave("../../figures/plots/aux_question_framing_comparison.pdf",
#        width = 6,
#        height = 6)

5.0.0.1.2 3D scatter plot of participant clusters
plot_ly(x = df.cluster_whs$whether,
        y = df.cluster_whs$how,
        z = df.cluster_whs$sufficient,
        type = "scatter3d",
        mode = "markers",
        color = as.factor(df.cluster_whs$cluster),
        colors = c("#e41a1c", "#377eb8")) %>% 
  layout(showlegend = FALSE,
         title = "Participant clusters",
         scene = list(
           xaxis = list(title = "whether"),
           yaxis = list(title = "how"),
           zaxis = list(title = "sufficient")))
5.0.0.1.3 Model fits for individual participants
df.plot = fit.brm_exp3_causal_whs %>% 
  fitted() %>% 
  as_tibble() %>% 
  select(prediction = Estimate) %>% 
  bind_cols(df.exp3.causal.long) %>% 
  relocate(prediction, .after = last_col())

df.text = df.plot %>% 
  group_by(participant) %>% 
  summarize(r = round(cor(prediction, rating), 2)) %>% 
  ungroup() %>% 
  mutate(r = str_c("r = ", r),
         prediction = 1,
         rating = 110)
`summarise()` ungrouping output (override with `.groups` argument)
ggplot(data = df.plot,
       mapping = aes(x = prediction,
                     y = rating)) + 
  geom_abline(intercept = 0, 
              slope = 1,
              color = "blue",
              alpha = 0.5) +
  geom_point(alpha = 0.3,
             size = 1) +
  geom_text(data = df.text,
            mapping = aes(label = r),
            hjust = 0,
            color = "red",
            size = 4) +
  facet_wrap(facets = vars(participant),
             ncol = 7) +
  labs(x = "model prediction",
       y = "causal responsibility rating") + 
  coord_cartesian(xlim = c(0, 100),
                  ylim = c(0, 100),
                  expand = F,
                  clip = "off") +
  scale_x_continuous(breaks = seq(0, 100, 25),
                     limits = c(0, 100)) +
  scale_y_continuous(breaks = seq(0, 100, 25),
                     expand = expansion(mult = c(0, 0))) +
  theme_classic() + 
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        text = element_text(size = 16),
        panel.spacing.x = unit(0.75, "cm"),
        panel.spacing.y = unit(1, "cm"),
        plot.margin = margin(t = 0.5, r = 0.5, unit = "cm"))

# ggsave("../../figures/plots/exp3_individual_scatter_plots.pdf",
#        width = 12,
#        height = 8)

5.0.0.1.4 Individual participant model fits
df.plot = df.exp3.causal.individual_fit %>% 
  select(participant, contains("r_"), contains("rmse_")) %>% 
  pivot_longer(cols = -participant,
               names_to = c("measure", "model"),
               names_sep = "_",
               values_to = "fit")

ggplot(data = df.plot %>% 
         filter(measure == "r"),
       mapping = aes(x = fit,
                     fill = model)) + 
  geom_density(alpha = 0.5) + 
  labs(x = "correlation value") + 
  scale_x_continuous(breaks = seq(0, 1, 0.2),
                     limits = c(0, 1)) + 
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) + 
  scale_fill_brewer(palette = "Set1") + 
  # ggplot2::theme_classic() + 
  theme(legend.position = c(0.3, 0.95),
        legend.direction = "horizontal",
        text = element_text(size = 20))

# ggsave("../../figures/plots/exp3_clusters_densities.pdf",
#        width = 8,
#        height = 6)

5.0.0.1.5 Clusters of participants
df.plot = df.exp3.causal.long %>% 
  left_join(df.cluster_whs %>% 
              select(participant, cluster) %>% 
              mutate(participant = as.numeric(participant),
                     cluster = as.factor(cluster)) %>% 
              group_by(cluster) %>% 
              mutate(label = factor(cluster, levels = 1:2, labels = str_c("n = ", group_size(.)))) %>% 
              ungroup(),
            by = "participant")

ggplot(data = df.plot,
       mapping = aes(x = ball,
                     y = rating,
                     group = label,
                     fill = label)) + 
  geom_point(shape = 21,
             alpha = 0.2,
             position = position_jitterdodge(jitter.width = 0.1,
                                             jitter.height = 0,
                                             dodge.width = 0.75)) + 
  stat_summary(fun.data = "mean_cl_boot",
               geom = "linerange", 
               mapping = aes(color = label),
               size = 1,
               position = position_dodge(width = 0.75)) + 
  stat_summary(fun = "mean",
               geom = "point", 
               size = 4,
               shape = 21,
               position = position_dodge(width = 0.75)) + 
  facet_wrap(facets = vars(clip), ncol = 8) + 
  labs(y = "causal responsibility rating") + 
  scale_fill_brewer(palette = "Set1") + 
  scale_color_brewer(palette = "Set1") + 
  theme(text = element_text(size = 24),
        legend.position = "bottom")

# ggsave("../../figures/plots/exp3_clusters_points.pdf",
#        width = 12,
#        height = 8)

5.0.0.1.6 Ternary plots for individual regression results
library("ggtern")
Registered S3 methods overwritten by 'ggtern':
  method           from   
  grid.draw.ggplot ggplot2
  plot.ggplot      ggplot2
  print.ggplot     ggplot2
--
Remember to cite, run citation(package = 'ggtern') for further info.
--

Attaching package: 'ggtern'
The following objects are masked from 'package:gridExtra':

    arrangeGrob, grid.arrange
The following objects are masked from 'package:ggplot2':

    aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob,
    ggsave, layer_data, theme_bw, theme_classic, theme_dark,
    theme_gray, theme_light, theme_linedraw, theme_minimal, theme_void
df.plot = df.exp3.causal.individual_fit %>% 
  select(participant, fit_whs) %>% 
  mutate(estimates = map(fit_whs, tidy)) %>% 
  select(participant, estimates) %>% 
  unnest(estimates) %>% 
  filter(str_detect(term, "b_"),
         !str_detect(term, "Intercept")) %>% 
  mutate(term = str_remove(term, "b_")) %>% 
  select(participant, term, estimate) %>% 
  pivot_wider(names_from = term,
              values_from = estimate) %>% 
  mutate_at(vars(whether, how, sufficient),
            .funs = list(norm = ~ . / (how + whether + sufficient))) %>% 
  mutate(color = 0,
         color = ifelse(how_norm == max(how_norm), 1, color),
         color = ifelse(whether_norm == max(whether_norm), 2, color),
         color = ifelse(sufficient_norm == max(sufficient_norm), 3, color),
         color = factor(color))

ggplot(data = df.plot,
       mapping = aes(x = how,
                     y = sufficient,
                     z = whether,
                     color = color)) + 
  geom_point(alpha = 0.7,
             size = 2,
             show.legend = F) + 
  scale_color_manual(values = c("black", "red", "green", "blue")) + 
  coord_tern() +
  theme_showarrows() + 
  theme(text = element_text(size = 20),
        tern.axis.title.T = element_blank(),
        tern.axis.title.L = element_blank(),
        tern.axis.title.R = element_blank())
Coordinate system already present. Adding new coordinate system, which will replace the existing one.
# ggsave(str_c("../../figures/plots/exp3_individual_regression_ternary_plot_scaled.pdf"),
#        width = 5,
#        height = 5)

6 Session info

sessionInfo()
R version 3.6.2 (2019-12-12)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggtern_3.3.0     forcats_0.5.0    stringr_1.4.0    dplyr_1.0.0     
 [5] purrr_0.3.4      readr_1.3.1      tidyr_1.1.0      tibble_3.0.1    
 [9] tidyverse_1.3.0  patchwork_1.0.0  egg_0.4.5        gridExtra_2.3   
[13] plotly_4.9.2.1   png_0.1-7        tidybayes_2.0.3  Hmisc_4.4-0     
[17] ggplot2_3.3.1    Formula_1.2-3    survival_3.1-12  lattice_0.20-41 
[21] janitor_2.0.1    jsonlite_1.6.1   xtable_1.8-4     corrr_0.4.2     
[25] brms_2.13.0      Rcpp_1.0.4.6     broom_0.5.6      lme4_1.1-23     
[29] Matrix_1.2-18    kableExtra_1.1.0 knitr_1.28      

loaded via a namespace (and not attached):
  [1] proto_1.0.0          tidyselect_1.1.0     htmlwidgets_1.5.1   
  [4] munsell_0.5.0        codetools_0.2-16     statmod_1.4.34      
  [7] DT_0.13              miniUI_0.1.1.1       withr_2.2.0         
 [10] Brobdingnag_1.2-6    colorspace_1.4-1     highr_0.8           
 [13] rstudioapi_0.11      stats4_3.6.2         robustbase_0.93-6   
 [16] bayesm_3.1-4         bayesplot_1.7.2      labeling_0.3        
 [19] emmeans_1.4.7        rstan_2.19.3         farver_2.0.3        
 [22] bridgesampling_1.0-0 coda_0.19-3          vctrs_0.3.1         
 [25] generics_0.0.2       TH.data_1.0-10       xfun_0.14           
 [28] R6_2.4.1             markdown_1.1         HDInterval_0.2.2    
 [31] assertthat_0.2.1     promises_1.1.0       scales_1.1.1        
 [34] multcomp_1.4-13      nnet_7.3-14          gtable_0.3.0        
 [37] processx_3.4.2       sandwich_2.5-1       rlang_0.4.6         
 [40] splines_3.6.2        lazyeval_0.2.2       acepack_1.4.1       
 [43] checkmate_2.0.0      inline_0.3.15        yaml_2.2.1          
 [46] reshape2_1.4.4       abind_1.4-5          modelr_0.1.8        
 [49] threejs_0.3.3        crosstalk_1.1.0.1    backports_1.1.7     
 [52] httpuv_1.5.3.1       rsconnect_0.8.16     tensorA_0.36.1      
 [55] tools_3.6.2          bookdown_0.19        ellipsis_0.3.1      
 [58] RColorBrewer_1.1-2   ggridges_0.5.2       latex2exp_0.4.0     
 [61] plyr_1.8.6           base64enc_0.1-3      ps_1.3.3            
 [64] prettyunits_1.1.1    rpart_4.1-15         zoo_1.8-8           
 [67] haven_2.3.1          cluster_2.1.0        fs_1.4.1            
 [70] magrittr_1.5         data.table_1.12.8    colourpicker_1.0    
 [73] reprex_0.3.0         mvtnorm_1.1-0        matrixStats_0.56.0  
 [76] hms_0.5.3            shinyjs_1.1          mime_0.9            
 [79] evaluate_0.14        arrayhelpers_1.1-0   shinystan_2.5.0     
 [82] jpeg_0.1-8.1         readxl_1.3.1         rstantools_2.0.0    
 [85] compiler_3.6.2       crayon_1.3.4         minqa_1.2.4         
 [88] StanHeaders_2.21.0-3 htmltools_0.4.0      mgcv_1.8-31         
 [91] later_1.1.0.1        lubridate_1.7.8      DBI_1.1.0           
 [94] dbplyr_1.4.4         MASS_7.3-51.6        boot_1.3-25         
 [97] compositions_1.40-5  cli_2.0.2            parallel_3.6.2      
[100] igraph_1.2.5         pkgconfig_2.0.3      foreign_0.8-76      
[103] xml2_1.3.2           svUnit_1.0.3         dygraphs_1.1.1.6    
[106] webshot_0.5.2        estimability_1.3     rvest_0.3.5         
[109] snakecase_0.11.0     callr_3.4.3          digest_0.6.25       
[112] rmarkdown_2.2        cellranger_1.1.0     htmlTable_1.13.3    
[115] shiny_1.4.0.2        gtools_3.8.2         nloptr_1.2.2.1      
[118] lifecycle_0.2.0      nlme_3.1-148         viridisLite_0.3.0   
[121] fansi_0.4.1          pillar_1.4.4         loo_2.2.0           
[124] fastmap_1.0.1        httr_1.4.1           DEoptimR_1.0-8      
[127] pkgbuild_1.0.8       glue_1.4.1           xts_0.12-0          
[130] shinythemes_1.1.2    stringi_1.4.6        blob_1.2.1          
[133] latticeExtra_0.6-29